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Abstract. We present numerical schemes for the strong solution of linear stochastic differen- 
tial equations driven by an arbitrary number of Wiener processes. These schemes are based on the 
Neumann (stochastic Taylor) and Magnus expansions. Firstly, we consider the case when the gov- 
erning linear diffusion vector fields commute with each other, but not with the linear drift vector 
field. We prove that numerical methods based on the Magnus expansion are more accurate in the 
mean-square sense than corresponding stochastic Taylor integration schemes. Secondly, we derive 
the maximal rate of convergence for arbitrary multi-dimensional stochastic integrals approximated 
by their conditional expectations. Consequently, for general nonlinear stochastic differential equa- 
tions with non-commuting vector fields, we deduce explicit formulae for the relation between error 
and computational costs for methods of arbitrary order. Thirdly, we consider the consequences in 
two numerical studies, one of which is an application arising in stochastic linear-quadratic optimal 
control. 
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1. Introduction. We are interested in designing efficient numerical schemes for 
the strong approximation of linear Stratonovich stochastic differential equations of 
the form 

yt = yo + Y. adr)yrdW^, (i.i) 

where y E MP, = t, {W^, . . . , W^) is a d-dimensional Wiener process and ao{t) 
and ttiit) are given p x p coefficient matrices. We call 'ao(i) y' the linear drift vector 
field and 'ai(t) y' for i — I, . . . , d the linear diffusion vector fields. We can express the 
stochastic differential equation ()l.ip more succinctly in the form 

y = yo + Koy, (1.2) 

where K = Kq + Ki + ■ • • + K^; and (K^ o y)t = /g ai(r) yr dW^. The solution of 
the integral equation for y is known as the Neumann series, Peano-Baker series, 
Feynman-Dyson path ordered exponential or Chen-Fleiss series 

yt = {\- K)-i o 2/0 = (I + K + K2 + + • • • ) o yo . 

The flow-map or fundamental solution matrix St maps the initial data yo to the 
solution yt = St yo at time t > 0. It satisfies an analogous matrix valued stochastic 
differential equation to ()1.2|) with the p x p identity matrix as initial data. The 
logarithm of the Neumann expansion for the flow-map is the Magnus expansion. We 
can thus write the solution to the stochastic differential equation (jl.ip in the form 

yt = (expert) yo, 
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where 

at = ln((l - K)-^ o /) = K o / + o / - i(K o /)2 + . . . . (1.3) 

See Magnus [41], Kunita [35], Azencott [3], Ben Arous [4|, Strichartz [54], Castell [12], 
Burrage [7], Burrage and Burrage [8] and Baudoin [6] for the derivation and conver- 
gence of the original and also stochastic Magnus expansion; Iserles, Munthe-Kaas, 
N0rsett and Zanna [3D] for a deterministic review; Lyons [35] and Sipilainen [5^ for 
extensions to rough signals; Lyons and Victoir [1^ for a recent application to prob- 
abilistic methods for solving partial differential equations; and Sussmann [5^ for a 
related product expansion. 

In the case when the coefficient matrices ai{t) — a.i, i — 0, . . . , d are constant and 
non-commutative, the solution to the linear problem (jl.ip is non-trivial and given by 
the Neumann series or stochastic Taylor expansion (see Kloeden and Platen [331) 

oo 

= X! X! Jat-a^ (t) aai--- VO , (1-4) 

where 

it) ^ f ■■■ /^"' dwi' ■ ■ ■ dwg^ dw^: ■ 

Jo Jo Jo 

Here is the set of all combinations of multi-indices a — {ai, . . . , ai} of length i 
with ak e {0, 1, . . . , 0?} for k — I, . . . , £. There are some special non-commutative cases 
when we can write down an explicit analytical solution. For example the stochastic 
differential equation dyt = aiyt <iWt + ytCi2 '^Wi with the identity matrix as initial 
data has the explicit analytical solution yt = exp(aiTy/) • exp(a2Wt'^). However in 
general we cannot express the Neumann solution series (|1.4p in such a closed form. 

Classical numerical schemes such as the Euler-Maruyama and Milstein methods 
correspond to truncating the stochastic Taylor expansion to generate global strong 
order 1/2 and order 1 schemes, respectively. Stochastic Runge-Kutta numerical meth- 
ods have also been derived — see Kloeden and Platen j33j and Talay |57j . At the linear 
level, the Neumann, stochastic Taylor and Runge-Kutta type methods are equivalent. 
In the stochastic context, Magnus integrators have been considered by Castell and 
Gaines Burrage [2|, Burrage and Burrage _8J and Misawa [46j . 

We present numerical schemes based on truncated Neumann and Magnus ex- 
pansions. Higher order multiple Stratonovich integrals are approximated across each 
time-step by their expectations conditioned on the increments of the Wiener processes 
on suitable subdivisions (see Newton [4^ and Gaines and Lyons [22]). What is new 
in this paper is that we: 

1. Prove the strong convergence of the truncated stochastic Magnus expansion 
for small stepsize; 

2. Derive uniformly accurate higher order stochastic integrators based on the 
Magnus expansion in the case of commuting linear diffusion vector fields; 

3. Prove the maximal rate of convergence for arbitrary multi-dimensional stochas- 
tic integrals approximated by their conditional expectations; 

4. Derive explicit formulae for the relation between error and computational 
costs for methods of arbitrary order in the case of general nonlinear, non- 
commuting governing vector fields. 
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Our results can be extended to nonlinear stochastic differential equations with anal- 
ogous conditions on the governing nonlinear vector fields, where the exponential Lie 
series (replacing the Magnus expansion) can be evaluated using the Castell-Gaines 
approach. 

In the first half of this paper, sections 2-5, we focus on proving the convergence of 
the truncated Magnus expansion and establishing Magnus integrators that are more 
accurate than Neumann (stochastic Taylor) schemes of the same order. The numerical 
schemes we present belong to the important class of asymptotically efficient schemes 
introduced by Newton [49'. Such schemes have the optimal minimum leading error 
coefficient among all schemes that depend on increments of the underlying Wiener 
process only. Castell and Gaines [131 E] prove that the order 1/2 Magnus integrator 
driven by a d-dimensional Wiener process and a modified order 1 Magnus integrator 
driven by a 1-dimensional Wiener process are asymptotically efficient. We extend 
this result of Castell and Gaines to an arbitrary number of driving Wiener processes. 
We prove that if we assume the linear diffusion vector fields commute, then an anal- 
ogously modified order 1 Magnus integrator and a new order 3/2 Magnus integrator 
are globally more accurate than their corresponding Neumann integrators. 

There are several potential sources of cost contributing to the overall computa- 
tional effort of a stochastic numerical integration scheme. The main ones are the 
efforts associated with: 

• Evaluation: computing (and combining) the individual terms and special 
functions such as the matrix exponential; 

• Quadrature: the accurate representation of multiple Stratonovich integrals. 
There are usually fewer terms in the Magnus expansion compared to the Neumann 
expansion to the same order, but there is the additional computational expense of 
computing the matrix exponential. When the cost of computing the matrix exponen- 
tial is not significant, due to their superior accuracy we expect Magnus integrators 
to be preferable to classical stochastic numerical integrators. This will be the case 
for systems that are small (see Moler and Van Loan [47] and Iserles and Zanna [3T] ) 
or for large systems when we only have to compute the exponential of a large sparse 
matrix times given vector data for which we can use Krylov subspace methods (see 
Moler and Van Loan [?7] and Sidje ^52]). Magnus integrators are also preferable when 
using higher order integrators (applied to non-sparse systems of any size) when high 
accuracies are required. This is because in this scenario, quadrature computational 
cost dominates integrator effort. 

In the second half of this paper, sections 6-8, we focus on the quadrature cost 
associated with approximating multiple Stratonovich integrals to a degree of accuracy 
commensurate with the order of the numerical method implemented. Our conclusions 
apply generally to the case of nonlinear, non-commuting governing vector fields. The 
governing set of vector fields and driving path process {W^, . . . , W^) generate the 
unique solution process y G MP to the stochastic differential equation Ijl.ip . For a 
scalar driving Wiener process W the Ito map W y is continuous in the topology of 
uniform convergence. For a d-dimensional driving processes with d > 2 the Universal 
Limit Theorem implies that the Ito map {W-^, . . . , W^) i-^ y is continuous in the p- 
variation topology, in particular for 2 < p < 3 (see Lyons |38 , Lyons and Qian [39] 
and Malliavin 43J). Since Wiener paths with d > 2 have finite p- variation for p > 2, 
approximations to y constructed using successively refined approximations to the 
driving path will only converge to the correct solution y if we include information 
about the Levy chordal areas of the driving path (the L^-norm of the 2-variation of 
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a Wiener process is finite though). Hence if we want to implement a scheme using 
adaptive stepsize we should consider order 1 or higher pathwise stochastic numerical 
methods (see Gaines and Lyons [22]). 

However simulating multiple Stratonovich integrals accurately is costly! For clas- 
sical accounts of this limitation on applying higher order pathwise stochastic numerical 
schemes see Kloeden and Platen ^31 p. 367], Milstein O p. 92] and Schurz [511 p. 58] 
and for more recent results see Gaines and Lyons [3TJ[22], Wiktorsson [55], Cruzeiro, 
MaUiavin and Thalmaier [T7], Stump and Hill [5S| and Giles [MIUS]. 

Taking a leaf from Gaines and Lyons |22j we consider whether it is computa- 
tionally cheaper to collect a set of sample data over a given time interval and then 
evaluate the solution (conditioned on that sample data), than it is to evaluate the 
solution frequently, say at every sample time. The resounding result here is that of 
Clark and Cameron [16] who prove that when the multiple Stratonovich integral J12 is 
approximated by its expectation conditioned on intervening sample points, the max- 
imal rate of L^-convergence is of order hlQ^I"^ where h is the integration steplength 
and Q is the sampling rate. We extend this result to multiple Stratonovich integrals 
Jax,...,ai of arbitrary order approximated by their expectation conditioned on inter- 
vening information sampled at the rate Q. Indeed we prove that the maximal rate 
of convergence is hf'l'^ jQ^I"^ when ai, . . . , are non-zero indices (and an improved 
rate of convergence if some of them are zero). In practice the key information is 
how the accuracy achieved scales with the effort required to produce it on the global 
interval of integration say [0, T] where T — Nh. We derive an explicit formula for 
the relation between the global error and the computational effort required to achieve 
it for a multiple Stratonovich integral of arbitrary order when the indices ai, . . . ,ai 
are distinct. This allows us to infer the effectiveness of strong methods of arbitrary 
order for systems with non-commuting vector fields. For a given computational effort 
which method delivers the best accuracy? The answer not only relies on methods 
that are more accurate at a given order. It also is influenced by three regimes for 
the stepsize that are distinguished as follows. In the first large stepsize regime the 
evaluation effort is greater than the quadrature effort; higher order methods produce 
superior performance for given effort. Quadrature effort exceeds evaluation effort in 
the second smaller stepsize regime. We show that in this regime when c? = 2, or when 
d > 3 and the order of the method M < 3/2, then the global error scales with the 
computational effort with an exponent of —1/2. Here more accurate higher order 
methods still produce superior performance for given effort; but not at an increas- 
ing rate as the stepsize is decreased. However when d > 3 for strong methods with 
M > 2 the global error verses computational effort exponent is worse than —1/2 and 
this distinguishes the third very small stepsize regime. The greater exponent means 
that eventually lower order methods will deliver greater accuracy for a given effort. 

We have chosen to approximate higher order integrals over a given time step by 
their expectations conditioned on the increments of the Wiener processes on suit- 
able subdivisions. This is important for adaptive time-step schemes (Gaines and 
Lyons [5^) and filtering problems where the driving processes (say and W"^) are 
observed signals. However it should be noted that Wiktorsson [58j has provided a 
practical method for efficiently sampling the set of multiple Stratonovich multiple in- 
tegrals {Jij'. i,j — 1, . . . , d} across a given time-step associated with a d-dimensional 
driving process (see Gilsing and Shardlow [23] for a practical implementation). Wik- 
torsson simulates the tail distribution in a truncated Karhunen-Loeve Fourier series 
approximation of these integrals which produces a convergence rate of order hd?/"^ /Q 
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where Q is analogously the number of required independent normally distributed 
samples. 

Other potential sources of computational effort might be path generation and 
memory access. Path generation effort depends on the application context. This cost 
is at worst proportional to the quadrature effort where we could subsume it. Memory 
access efforts depend on the processing and access memory environment. To reveal 
higher order methods (which typically require more path information) in the best 
light possible, we have ignored this effect. 

Our paper is outlined as follows. We start in Section 2 by proving that the 
exponential of every truncation of the Magnus series converges to the solution of our 
linear stochastic differential equation Ijl.ip . In Section 3 we define the strong error 
measures we use and how to compute them. Using these, we explicitly compare the 
local and then global errors for the Magnus and Neumann integrators in Section 4 
and thus establish our stated results for uniformly accurate Magnus integrators. In 
Section 5 we show that when the linear diffusion vector fields do not commute we 
cannot expect the corresponding order 1 Magnus integrator to in general be globally 
more accurate than the order 1 Neumann integrator. We then turn our attention in 
Section 6 to the method of approximating multiple Stratonovich integrals by their 
conditional expectations. We prove the maximal rate of convergence for an arbitrary 
multiple Stratonovich integral in Section 6. We then use this result in Section 7 to 
show how the global error scales with the computational effort for numerical schemes 
of arbitrary order. The shuffle algebra of multiple Stratonovich integrals generated 
by integration by parts allows for different representions and therefore bases for the 
solution of a stochastic differential equation. Some choices of basis representation are 
more efficiently approximated than others and we investigate in Section 8 the impact of 
this choice. In Section 9 we present numerical experiments that reffect our theoretical 
results. To illustrate the superior accuracy of the uniformly accurate Magnus methods 
we apply them to a stochastic Riccati differential system that can be reformulated as a 
linear system which has commuting diffusion vector fields. Since for the linear system 
expensive matrix-matrix multiplications can be achieved independent of the path, the 
Neumann method performs better than an explicit Runge-Kutta type method applied 
directly to the nonlinear Riccati system. We also numerically solve an explicit linear 
system with governing linear vector fields that do not commute for two and also 
three driving Wiener processes — Magnus integrators also exhibit superior accuracy 
in practice in these cases also. Lastly in Section 10, we outline how to extend our 
results to nonlinear stochastic differential equations and propose further extensions 
and applications. 

2. Strong convergence of truncated Magnus series. We consider here the 
case when the stochastic differential equation (jl.ip is driven by d Wiener processes 
with constant coefficient matrices ai(t) = Oi, i = 0, 1, . . . , d. The Neumann expansion 
has the form shown in p.4p . We construct the Magnus expansion by taking the 
logarithm of this Neumann series as in (jl.3p . In Appendix \X\ we explicitly give the 
Neumann and Magnus expansions up to terms with L^-norm of order 3/2. Let at 
denote the truncated Magnus series 



where Q„i denotes the finite set of multi- indices a for which || Jq||l2 is of order up to 
and including i™. Note that here m is a half-integer index, m = 1/2, 1, 3/2, . . .. The 




(2.1) 
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terms are linear combinations of finitely many (more precisely exactly length a) 
products of the ai, i = 0,1, . . . ,d. Let |Qm| denote the cardinality of Q^. 

Theorem 2.1 (Convergence). For any t < 1, the exponential of the truncated 
Magnus series, expat, is square-integrable. Further, ifyt is the solution of the stochas- 
tic differential equation there exists a constant C{m) such that 

||2/t - expat • 2/0IL2 <C^(m)t™+^/^ (2.2) 



Proof. First we show that expat e Li^ . Using the expression (|2.ip for at, we see 
that for any number k, {(Tt)^ is a sum of |Qm|'^ terms, each of which is a fc-multiple 
product of terms JaCa- It follows that 

||('7t)''||^2 < ( l|Ca||op) • ^ \\Ja{l)Ja{2) 

Z = l,...,fe 

Note that the maximum of the operator norm || • ||op of the coefficient matrices 
is taken over a finite set. Repeated application of the product rule reveals that 
the product Ja{i)Ja{i)^ where a(i) and a{j) are multi-indices of length t(a(i)) and 
(.{a{j)), is a linear combination of 2^'^"(*^)+^("'^^^)~^ multiple Stratonovich integrals. 
Since i[a{i)) < 2m for i = I, . . . , k, each term ' Jq(i) Jq(2) • • • Ja(k) in (|2.3p is thus the 
sum of at most 2^™'^^^ Stratonovich integrals J/3. We also note that k < £{(]) < 2mk. 

From equation (5.2.34) in Kloeden and Platen [33], every multiple Stratonovich 
integral J/3 can be expressed as a finite sum of at most 2^^'^^^^ multiple Ito integrals 
Ij with < e{P). Further, from Remark 5.2.8 in Kloeden and Platen [33], £(7) + 
11(7) > + n(/9), where n{(3) and n(7) denote the number of zeros in (3 and 7, 
respectively. From Lemma 5.7.3 in Kloeden and Platen [33], 

Noting that < e{l3) < 2mk and £(7) + n(7) > k, it follows that for t < 1, we 
have ||J/3||l2 < 2'*'"'^"^ i'^/^. Since the right hand side of equation (|2.3p consists of 
IQml'^ 2^™''^-'^ Stratonovich integrals Jp, we conclude that, 




Hence exp at is square-integrable. 

Second we prove (|2.2p . Let yt denote Neumann series solution (|1.4p truncated to 
included terms of order up to and including i™. We have 

\\yt - expat • yo||i2 < \\yt - yt\\^2 + \\yt - expat • 2/o||i2 • (2-4) 

We know yt G (see Gihman and Skorohod [53] or Arnold ,2 ). Furthermore, for 
any order m, yt corresponds to the truncated Taylor expansion involving terms of 
order up to and including i™. Hence yt is a strong approximation to yt to that order 
with the remainder consisting of C'(i™+^/^) terms (see Proposition 5.9.1 in Kloeden 
and Platen 33J). It follows from the definition of the Magnus series as the logarithm 
of the flow-map Neumann series, that the terms of order up to and including in 
expat • yo correspond with yt; the error consists of 0(i™+^/^) terms. □ 



•••Ja(fc)llL- (2.3) 
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Convergence of approximations based on truncations of the stochastic Taylor 
expansion has been studied in Kloeden and Platen see Propositions 5.10.1, 5.10.2, 
and 10.6.3. Ben Arous [1| and Castell [T^ prove the remainder of the exponential of 
any truncation of the Magnus series is bounded in probability as t — s- (in the full 
nonlinear case). Burrage [7 shows that the first terms up to and including order 3/2 
Magnus expansion coincide with the terms in the Taylor expansion of the same order. 
Our result holds for any order in for sufficiently small t. A more detailed analysis 
is needed to establish results concerning the convergence radius. Similar arguments 
can be used to study the non-constant coefficient case with suitable conditions on 
the coefficient matrices (see Proposition 5.10.1 in Kloeden and Platen [33] for the 
corresponding result for the Taylor expansion). 

Note that above and in subsequent sections, one may equally consider a stochastic 
differential equation starting at time > with square- integrable J-'tj, -measurable 
initial data t/o- Here {J-'t)t>o denotes the underlying filtration. 

3. Global and local error. Suppose S't„,t„+i and S't„,t„+i are the exact and 
approximate flow-maps across the interval [i„,i„_|_i], respectively; both satisfying the 
usual flow-map semi-group property: composition of flow- maps across successive inter- 
vals generates the flow-map across the union of those intervals. We call the difference 
between the exact and approximate flow-maps 



Rt 



(3.1) 



the local flow remainder. For an approximation yt„^-i across the interval [i„, the 
local remainder is thus Rt„,tn+iyt„- Our goal here is to see how the leading order 
terms in the local remainders accumulate, contributing to the global error. 

Definition 3.1 (Strong global error). We define the strong global error associ- 
ated with an approximate solution yx to the stochastic differential equation (|l.ip over 
the global interval of integration [0,T] as £ = \\yT — yxWh"^- 

The global error can be decomposed additively into two components, the global 
truncation error due to truncation of higher order terms, and the global quadrature 
error due to the approximation of multiple Stratonovich integrals retained in the 
approximation. If [0, T] 
h^T/N we have 



^n=o\^n,tn+i\ whcre tn = nk then for small stepsize 



£ = 



( n n 

\ri=W-l n=N-l / 

I X! ^t„+i,t,^Rt„.t„+iSto,tr, I yo 



yo 



L2 



+ 0{max\\Rt„,t^^,\\%^h-^/^) . (3.2) 



The local flow remainder has the following form in the case of constant coefficients 
ai, i = 1, . . . , d, (see for example the integrators in Appendix [X]) : 



Rt 



^1 — ^ ^ J a {tn 1 ) Cq 



Here a is a multi-index and the terms Cq represent products or commutations of the 
constant matrices a^. The Jq represent Stratonovich integrals (or linear combinations, 
of the same order, of products of integrals including permutations of a). The global 
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error at leading order in the stepsize is thus 

n+l;tN 

Hence in the global truncation error we distinguish between the diagonal sum con- 
sisting of the the first sum on the right-hand side above, and the off-diagonal sum 
consisting of the second sum above with n ^ m. 

Suppose we include in our integrator all terms with local i^-norm up to and 
including 0{h^^). The leading terms Jq,Cq in Rt„.t„+i thus have i^-norm C(/i^^+^/^). 
Those with zero expectation will contribute to the diagonal sum, generating 0{h^^) 
terms in the global error, consistent with a global order M integrator. However those 
with with non-zero expectation contribute to the off-diagonal double sum. They will 
generate 0{h^'^~'^/'^) terms in the global error. We must thus either include them in 
the integrator, or more cheaply, only include their expectations — the corresponding 
terms (Jq, — E(Jq))cq. of order in Rt„.t„+i will then have zero expectation 

and only contribute through the diagonal sum — see for example Milstein [44|, p. 12]. 
This also guarantees the next order term in the global error estimate p.2p . whose 
largest term has the upper bound max„ ||-Rt„,t„+i ||^2^ h^^^'^, only involves higher order 
contributions to the leading 0{h'^^) error. 

Note that high order integrators may include multiple Stratonovich integral terms. 
We approximate these multiple integrals by their conditional expectations to the local 
order of approximation /i*^"*"^/^ of the numerical method. Hence terms in the inte- 
grator of the form JaCa are in fact approximated by K{Ja\J^Q)ca, their expectation 
conditioned on intervening path information !Fq (see Section |6] for more details). This 
generates terms of the form (Jq — E(Jq|JFq))cq in the local flow remainder, which 
have zero expectation and hence contribute to the global error through the diagonal 
sum generating 0{h^) terms. 

4. Uniformly accurate Magnus integrators. Our goal is to identify a class 
Magnus integrators that are more accurate than Neumann (stochastic Taylor) integra- 
tors of the same order for any governing set of linear vector fields (for the integrators 
of order 1 and 3/2 we assume the diffusion vector fields commute). We thus com- 
pare the local accuracy of the Neumann and Magnus integrators through the leading 
terms of their remainders. We consider the case of constant coefficient matrices a^, 
i = 0,l,...,d. 

The local fiow remainder of a Neumann integrator i?"*^" is simply given by the 
terms not included in the flow-map Neumann approximation. Suppose a is the trun- 
cated Magnus expansion and that p is the corresponding remainder, i.e. a — a -\- p. 
Then the local flow remainder R'^'^^ associated with the Magnus approximation is 

i^'n'^s = exp cr - exp a 

= exp ((7 -\- p) — exp & 

= p+^{ap + pa) + o{a'p). (4.i) 

Hence the local flow remainder of a Magnus integrator i?™'*^ jg the truncated Magnus 
expansion remainder p, and higher order terms ^ {ap + pa) that can contribute to the 
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global error at leading order through their expectations. For the integrators considered 
in this section these higher order terms do not contribute in this way, however for the 
order 1 integrator we consider in the next section they do. 

Definition 4.1 (Uniformly accurate Magnus integrators). When the linear dif- 
fusion vector fields commute so that [ttj, aj] = for all i, j ^ 0, we define the order 1 
and order 3/2 uniformly accurate Magnus integrators by 

d 

i=l 

and 

d 

^tlitl+i = Jocio + ^{JiCii + ^{Jio - Joi)[ao,ai\ + ^[di, K,oo]]) • 

i=l 

By uniformly we mean for any given set of governing linear vector fields (or 
equivalently coefficient matrices a^, i = 0,1,..., d) for which the diffusion vector 
fields commute, and for any initial data yo. 

Theorem 4.2 (Global error comparison). For any initial condition yo and suffi- 
ciently small fixed stepsize h — tn+i — tn, the order 1/2 Magnus integrator is globally 
more accurate in than the order 1 /2 Neumann integrator. If in addition we assume 
the linear diffusion vector fields commute so that [ai,aj] = for all i, j ^ 0, then the 
order 1 and 3/2 uniformly accurate Magnus integrators are globally more accurate 
in LP' than the corresponding Neumann integrators. In other words, if denotes 
the global error of the order 1/2 Magnus integrator or the uniformly accurate Magnus 
integrators of order 1 or order 3/2, respectively, and f is the global error of the 
Neumann integrators of the corresponding order, then at each of those orders, 

^mag < ^neu _ (4 2) 

Proof. Let i?™'*^ ^nd R^^^ denote the local flow remainders corresponding to the 
Magnus and Neumann approximations across the interval [tni^n+i] with tn = nh. A 
direct calculation reveals that 

]E((^neu^T^„eu^ = E((i?'"^S)T^mag^ + Dm + 0{h^^+^''^) , 

where if we set R = i?""^" - iJ^^s then 

Dm = E{FFR"'''^) + E((i?">^s)T^^ ^ E{R^R) . (4.3) 

We now show explicitly in each of the three cases of the theorem that at leading order: 
^mag g^j^^ ^ uncorrelated and hence Dm is positive semi-definite. This implies 
that the local remainder for the Neumann expansion is larger than that of the Magnus 
expansion. Hereafter assume the indices i,j,k,l € {1,. . . ,d}. 

For the order 1/2 integrators we have to leading order (see the form of the ex- 
pansions in Appendix A) 

R ^ = 2 ("^y ~ ^i] ) 

i<j 
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and 



which arc uncorrclatcd by direct inspection. 

We henceforth assume [ai,aj] = for all For the uniformly accurate order 1 
integrator we have to leading order (again see Appendix A) 

i 

and 

R = i/i^tto + "^^{Jio + Joi){aoai + aiflo) + 3/1^(0? ao + aoai) 

i 

which again, by direct inspection are uncorrclatcd. 

For the uniformly accurate order 3/2 integrator the local flow remainders are 

^neu _ y~^(j^^Q _ jh^^aoaf + Jioittiaoai + {Jou — \h^)a^ao 

i 

+ y~^f (Joji + Joij)(li(ljO'0 + JjOiO'iO'OO'j 
i<j 

+ JiOjO'jC'OO'i + {JjiO + Jijo)0'00'iC'j^ 

+ ^ {Jijki -^{Jijki))aiakajai , 



ijykyl 



and 



^mag = ^ ^( J2 Jo - /l2 - 6 JiOi)[ai, [fli, Oq]] 

+ ^H'^oJiJj -HJiOj + Jjo^))[aj,[ai,ao]]■ 
i<j 

Consequently we have 

R = 3^ ^((6JioJi - JfJo - 2h'^)aoa1 + 2{jfJo - h^)aiaoai 

i 

+ {QJoiJi - jfJo - 2/i2)a2ao) 
+ g^^l {^{JiJoj + JjJoi) — JoJiJj)aiajao 

i<j ^ 

+ {JoJiJj + 3(JiOj" — JjOi))o,jO.OO.i 

+ [JaJiJj + 3{Jjoi — Jioj))aiaf)aj 

+ (3(JjJjO + JjJio) — JoJiJj)ciO<li<lj 

+ ^ {Jijki - E{Jijki))aiakajai . 

i,j,k,l 
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First we note that the terms in R of the form J^jkiaiakajai, for which at least three 
of the indices distinct, are uncorrelated with any terms in R™^^. We thus focus on 
the terms of this form with at most two distinct indices, namely 

i<j i=ij i 

and other remaining terms in R. Since E[Jioi| Ji] = h{jf — h)/6 and E[Jioj|Ji, Jj] = 
hJiJj/6 for i j the following conditional expectations are immediate 

E[(j2 Jo - h^){jfJo ~h^- 6J^0i)\J^] = , 

E[(J..,. - lh^)iJfJo -h^- 6J^0^)\J^] = , 

E[(J„; J,, - i/l')(jf Jo -h^~ GJ^O^)\J^, Jj] = , 

E[( Jo Ji Jj + 3( JiOj" — JjQi)) {JoJiJj ~ 3( JiOj + JjOi)) \ Ji^ — ^ ; 

E[( Jii Jjj — jh^){JoJiJj — 3( JiOj + JjOi))\Ji, Jj] = , 

E[Jiii Jj [JoJiJj — 3{JiOj + JjOi)) \ Ji, Jj] — . 

Hence the expectations of the terms shown are also zero. Secondly, direct computation 
of the following expectations reveals 

E((6Jo. J. - Jf Jo - 2h^)iJ^Jo -h^- 6J,o^)) = , 

E^(3(Ji Joj + JjJoi) — JoJiJj) {JoJiJj — 3( JiOj" + J?oi))) — . 

Hence R™'^^ and R are uncorrelated. 

The corresponding global comparison results (|4.2|) now follow directly in each of 
the three cases above using that the local errors accumulate and contribute to the 
global error as diagonal terms in the standard manner described in detail at the end 
of Section [3] Note that we include the terms j^h^lui, [a^, ao]] in the order 1 uniformly 
accurate Magnus integrator. These terms appear at leading order in the global re- 
mainder where they would otherwise generate a non-positive definite contribution to 
Dm in gS]). □ 

Note that the Magnus integrator ct'"'^' is an order 1 integrator without the terms 
j^h'^[ai,[ai,ao]]. However they are cheap to compute and including them in the 
integrator guarantees global superior accuracy independent of the set of governing 
coefficient matrices. 

5. Non-commuting linear diffusion vector fields. What happens when the 
linear diffusion vector fields do not commute, i.e. we have [ai,aj] 7^ for non-zero 
indices i ^ p. Hereafter assume A: e Consider the case of order 1 

integrators. The local flow remainders are 

ofloai + Joifliao) + ^ Jk^iaia^ak 

and 

^mag ^ ^ 1 (j^^ _ Joi)[aO, a.i] + ^{J^^J - \JiJ,j + ^Jpj) [Oj, flj]] 
i i^j 

+ ^ {{Jijk + \JjJki + \JkJij ^ ^JiJjJk)[o,iT[aj,ak]] 
i<j<k 

+ [Jjik + \ JiJkj + \ JkJji — \ JiJjJk)[aj, [a.j,afc]]) . 
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Computing Dm in (14. 3p gives 

Dm = h'Y. UujBUuj + h' U^kCU^jk 

where Uuj = (ajOi, aiajUi, afaj, Oj, aoaj, ajao)"^ G M^^^p and in addition we have that 
Uijk = {akO'jai,ajakai,akaiaj,aiakaj,ajaiak,aiajak)'^ G M^p^p. Here B,C £ M^pxSp 
consist of p X p diagonal blocks of the form bkelp and Ckelp where 



and 
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18 
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10 
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10 
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10 


31 


18 


24 


12 


144 


18 





18 


60 


36 


36 






12 





24 


36 


36 


36 






i^24 





12 


36 


36 


36/ 




/4 


1 


1 


1 


1 


-2 




1 


4 


1 


-2 


1 


1 




1 


1 


4 


1 


-2 


1 


6 


1 


-2 


1 


4 


1 


1 




1 


1 


-2 


1 


4 


1 






-2 


1 


1 


1 


1 


4 



Again, c is positive semi-definite with eigenvalues {^,^,^,^,0,0}. However b has 
eigenvalues {^,^^(5 + ^41), (5 - ^41), 0.94465, 0.0943205, -0.03897}, where the 
final three values are approximations to the roots of 288a;'^ — 2888a:^ + 14a; + 1. 

The eigenvalues of b and c, respectively, are multiple eigenvalues for the matrices 
B and C, respectively. This implies that there are certain matrix combinations and 
initial conditions, for which the order 1 Taylor approximation is more accurate in the 
mean-square sense than the Magnus approximation. However, the two negative values 
are small in absolute value compared to the positive eigenvalues. For the majority 
of systems, one can thus expect the Magnus aproximation to be more accurate (as 
already observed by Sipilainen [33], Burrage [7] and Burrage and Burrage [S]). For 
any given linear system of stochastic differential equations, the scheme that is more 
accurate can be identified using the results above. 

In this case there are terms from ^{ap+pa) in Magnus remainder expression (|4.ip 
that appear at leading order in the local flow remainder. These terms are of the form 
j^(ai[aj, [aj,ai]] -I- aj[ai, [ai,aj]])h'^ . They make a negative definite contribution to 
global error, though this can be negated by including them as cheaply computable 
terms in the Magnus integrator (indeed we recommend doing so). 

6. Quadrature. We start by emphasizing that there are two inherent scales: 

• Quadrature scale At on which the discrete Wiener paths are generated; 

• Evaluation scale h on which the stochastic differential equation is advanced. 
The idea is to approximate multiple Stratonovich integrals by their corresponding 
expectations conditioned on the cr-algebra representing intervening knowledge of the 
Wiener paths (Clark and Cameron [T^]; Newton [13]; Gaines and Lyons [22 ). Hence 
we approximate J„(i„,t„+i) by E ( Ja(i„, t„+i)| .fg) where 

Tq ^ {AWl^^^,: I = 1, . . . ,d; q ^ 0, . . . ,Q ^ I; n - 0, . . . , TV - 1} , 
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with AVFj*^ = _|_(q+i)At ~ ^t*„+gAt' = Q ^^e number of 

Wiener increments. We extend the result of Clark and Cameron 1 16J on the maximum 
rate of convergence to arbitrary order multiple Stratonovich integrals. 

Lemma 6.1 (Quadrature error). Suppose at least two of the indices in the multi- 
index a = {«!, • ■ • , ae} are distinct. Define 

j* = mm{j < e - 1: at ^ ai ,\fi,l > j} . 

Let Ti{a) denote the number of zeros in a, and n*(Q;) — n({Q!j. , . . . ,ai}). Then the 
error in the multiple Stratonovich integral approximation E (^Ja{tn,tn+i) \ J^q) is 

/ f^{e+n{a))/2 \ 

\\Ja{tn,tn+i) - E { Ja{tn, tn+l)\ ^q) 11^2 = ^ ( Q(n' (a) + l)/2 ) ' ^^■^'^ 



Proof. For any s,t, we write for brevity 

Ja{s,t)^E{j^{s,t)\J^Q) . 

We define as the multi-index obtained by deleting the last k indices, that is 
a~'' = {ai, • • • , ag-k}. We set = t„ + qAt, where g = 0, . . . , Q — 1. Then 



Ja {tjl ; t' 



n+l) 



q=Q " ^1 

Q-i n-i \ 

q=0 \k=l J 



Thus we have 



Ja{tji^tji^l^ J a{tn ^t^i-if-x) 

= ^ ^ I ^ ] {Ja-'' {tn, Tq) — J^^-k {t^, Tq)^ Jai,_f._i_j^....^ai ijq, T^+l) 
g=0 \fc=l 

+ ^Q-fc {tm Tq) (^Jai_k+i,...,at (''"?J ''"g+l) ~ Ja(_k+l,---,at i'^q^ '''9+1)) 
+ Ja{Tq, Tq+l) - Jaijq, Tg+l) . 



(6.2) 



We prove the assertion by induction over For 1 = 2^ the first two terms in the 
sum (j6.2p are zero and 



'Y{ja{Tq,Tq+l) - Ja{Tq,Tq+l)) 
9=0 



o(g(At)^+"(") 

Qe+n{a)-l 



= c 

Here we have used that fact that (. — 2, j* = 1 and thus ii{a) ~ n*(a). 
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Assume now that £ > 2. We will investigate the order of each of the three types of 
terms in (|6.2p separately. If at least two indices in a~'' are distinct, then by induction 
hypothesis 



12 \^ ^n*(Q-'=) + l 



and since — i„ = qAt and Q At = h, we have for each fc = l,...,^— 1, 

^ ^ II { Ja-'' {tn, Tq) — J^-k [tn, Tg)) Jae-k+l,---,at) {Tq^ ''"ij+l) II2 



Note that we have 

n*{a-'') + n{a) - n(a"''') + k>n*{a) + k> n*{a) + 1 . (6.3) 

If all indices in a^'^ are equal, then || Jq-i- (tn, r^) — •/q-': (^n, Tg)!^ = 0. 
For the second term in (16.21) we have for each k = 1 i — 1, 



^ - 2 

^ ' II Jq-S: (^n, Tq) { Jat-k + i,---,ae ('''g; ''"g+l) ^ "^af-fc + i,...,Q!f (''"g? '''g+l)) II2 
g=0 

, if . . . = ai-k+1 , 

O Q„(„)-l|r-°'=')+.-i ) ' otherwise . 

Note that in the second of these cases 

n(a) -n(a"'^') + fc- 1 > n*(a) + 1. (6.4) 
For the third term we have 

2^ Pa(T,,T,+ i) - Ja(r,,Tg+i)||2 = O ( 
g=0 



n(Q)-l 



Again, note that we have 

^ + n(a) - 1 > n*(a) + 1. (6.5) 



Equality holds in at least one of ()6.3|) to ()6.5p . To see this we distinguish the 
case when the last two indices are equal, a^-i = a^, and the case when the last two 
indices are distinct, ai-i ^ a(. If ai-i = ai, then 

n*(a-i) + n(a) - n(a-i) = n*(a) . 

Since in this case at least two indices in a"^ are distinct, equality holds for fc = 1 
in (|6.3p . If a^-i then 

n*(a) = n(Q!) — n(Q;~^) , 
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and thus equality holds for fc = 2 in (|6.4p . Hence the lemma follows. □ 
Note that each multiple Stratonovich integral Jaitmtn+i) can be thought of as 
an ^-dimensional topologically conical volume in {W^, . . . , Vl^^)-space. The surface of 
the conical volume is panelled with each panel distinguished by a double Stratonovich 
integral term involving two consecutive indices from a. The edges between the panels 
are distinguished by a triple integral and so forth. The conditional expectation ap- 
proximation K(^Ja{tn,tn+i)\ J'q) Can also be decomposed in this way. In the L^-error 
estimate for this approximation, the leading terms are given by sums over the panels 
which also confirm the estimate (I6.ip . 

Approximations of multiple Stratonovich integrals constructed using their condi- 
tional expectations are intimately linked to those constructed using paths Wl that 
are approximated by piecewise linear interpolations of the intervening sample points. 
The difference between the two approaches are asymptotically smaller terms. For 
more details see Wong and Zakai [53] , Kloeden and Platen [33 , Hofmann and Miiller- 
Gronbach f29l and Gyongy and Michaletzky [27] . 

7. Global error vs computational effort. We examine in detail the step- 
size/accuracy regimes for which higher order stochastic integrators are feasible and 
also when they become less efficient than lower order schemes. In any strong simula- 
tion there are two principle sources of computational effort. Firstly there is evaluation 
effort W™' associated with evaluating the vector fields, their compositions and any 
functions such as the matrix exponential. Secondly there is the quadrature effort U'^^^'^ 
associated with approximating multiple stochastic integrals to an accuracy commen- 
surate with the order of the method. For a numerical approximation of order M the 
computational evaluation effort measured in flops over N — Th~^ evaluation steps is 

U'^^' = {cMP^ + CE)Th~K 

Here p is the size of the system, cm represents the number of scalar-matrix multipli- 
cations and matrix-matrix additions for the order M truncated Magnus expansion, 
and cb is the effort required to compute the matrix exponential. Note that if we im- 
plement an order 1/2 method there is no quadrature effort. Hence since £ — 0{h^^^) 
we have £ = C'((W°™')"^^^)- 

Suppose we are required to simulate Jai - a^ (^n, ^n+i) with all the indices distinct 
with a global error of order h^^ ; naturally £ > 2 and M > £/2. 

Lemma 7.1. The quadrature effort hi measured in flops required to approximate 
Jai---ae(tn, tn+i) with a global error of order h^^ when all the indices are distinct and 
non-zero, is to leading order in h: 

where 

(3{M,£) ^ {£ - 1){2M + I - £) + 1 . 

Since we stipulate the global error associated with the multiple integral approximation 
to be £ = 0{h^'^), we have 

£ = o{U-^"^^^'^'^) . 

Proof. The quadrature effort required to construct E( jQj...Qj,(i„, J^g), 
which is a (£ - l)-multiple sum, over [0, T]isU^ 0{Q'^-^N) with = Th-^. Using 
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Table 7.1 

Slopes of the logarithm of the global error E verses the logarithm of the quadrature effort U, 
i.e. the exponent —M/l3{M,£), for different values of £ and M when Jai-.-af has distinct non-zero 
indices. 





e = 2 


e = 3 


£ = 4: 


£=5 £=e 


no zero index 




M = 1 


-1/2 








M = 3/2 


-1/2 


-1/2 






M = 2 


-1/2 


-2/5 


-1/2 




M = 5/2 


-1/2 


-5/14 


-5/14 


-1/2 ... 


M = 3 


-1/2 


-1/3 


-3/10 


-1/3 -1/2 



Lemma 16.11 to achieve a global accuracy of order ft,*^ and therefore local i^-norm of 
order for this integral, requires that Q — h^^^~^^^ . □ 

In Table 7.1 we quote values for the exponent —M/f3{M,£) for different values of 
M and £ in the case when all the distinct indices ai, . . . ,ae are non-zero. 

Suppose we are given a stochastic differential equation driven by a d-dimensional 
Wiener process with non-commuting governing vector fields. To successfully imple- 
ment a strong numerical method of order M we must guarantee that the global error 
associated with each multiple integral present in the integrator that is approximated 
by its conditional expectation is also of order M. If we implement a numerical method 
of order M < d/2, we will in general be required to simulate multiple Stratonovich 
integrals with distinct indices of length £ with 2 < £ < 2M < d. We will also have to 
simulate multiple integrals with repeated indices of length £r < 2M. These integrals 
will require the same or fewer quadrature points than those with distinct indices as we 
can take advantage of the repeated indices — see Section [8] for more details. Such in- 
tegrals therefore represent lower order corrections to the quadrature effort. Similarly 
multiple integrals with distinct indices that involve a zero index have index length 
that is one less than similar order multiple integrals with distinct non-zero indices. 
Hence they will also represent lower order corrections to the quadrature effort. 

To examine the scaling exponent in the relation between the global error £ and the 
quadrature effort VC^^^, which is the sum total of the efforts required to approximate 
all the required multiple integrals to order , we use Table 7.1 as a guide for the 
dominant scalings. For methods of order M < d/2, ii d = 2 and we implement a 
method of order M — 1, then the dominant exponent is —1/2. Similarly if d = 3 then 
order 1 and 3/2 methods also invoke a dominant scaling exponent of —1/2 for the 
integrals of length £ = 2 and £ = 3. If d = 4 then methods of order 1 and 3/2 have 
the same scaling exponent of —1/2, however the method of order 2 involves multiple 
integrals with three indices which are all distinct and the dominant scaling exponent 
for them is —2/5. 

If we implement a method of order M > d/2 then we will be required to simulate 
multiple Stratonovich integrals with distinct indices of length £ with 2 < £ < d. We 
must also simulate higher order multiple integrals with indices of length £r involving 
repetitions with d < £r < 2M; these may be cheaper to simulate than multiple 
integrals of the same length with distinct indices (again see Section [8|) . When d = 2 
the dominant scaling exponent is —1/2 for all orders. For c? > 3 the dominant scaling 
exponent is at best —1/2, and so forth. 

Lastly, we give an estimate for the critical stepsize ha- below which the quadrature 
effort dominates the evaluation effort. Since £ = M + 1 minimizes f3{M,£) we have 



EfRcient stochastic integrators 



17 



the following estimate. 

Corollary 7.2. For the case of general non-commuting governing vector fields 
and a numerical approximation of order AI , we have > ^i"'^'^ if and only if 

h > her where the critical stepsize 

/icr = 0\{T{cmP +Ce)) 1 , 

where ^max = maxjii, Af + 1}. 

In practice when we implement numerical methods for stochastic differential equa- 
tions driven by a c?-dimensional Wiener process we expect that for h > h^ the evalu- 
ation effort dominates the compuational cost. In this scenario integrators of order M 
scale like their deterministic counterparts. Consider what we might expect to see in a 
log-log plot of global error verses computational cost. As a function of increasing com- 
putational cost we expect the global error for each method to fan out with slope — Af , 
with higher order methods providing superior accuracy for a given effort. However 
once the quadrature effort starts to dominate, the scaling exponents described above 
take over. When d — 2 for example and all methods dress themselves with the scaling 
exponent —1/2, then we expect to see parallel graphs with higher order methods still 
providing superior accuracy for a given cost. However higher order methods that 
assume a scaling exponent worse than —1/2 will eventually re-intersect the graphs of 
their lower order counterparts and past that regime should not be used. 

Note that in the case when all the diffusion vector fields commute, methods of 
order 1 do not involve any quadrature effort and hence £ — ©((Z^"™')^^) . Using 
Lemma 16.11 we can by analogy with the arguments in the proof of Lemma 17.11 de- 
termine the dominant scaling exponents for methods of order M > 3/2. For example 
the L^-error associated with approximating Jgi by its expectation conditioned on in- 
tervening information is of order h^/'^/Q. Hence we need only choose Q — to 
achieve to achieve a global error of order 3/2. In this case the dominant scaling expo- 
nent is —1. However the L^-error associated with approximating Joij for i ^ j is of 
order h^/Q^^^ whereas for JiQj and J^o it is of order h^/Q. For the case of diffusing 
vector fields we do not need to simulate Joij, and so for a method of order 2 the 
dominant scaling exponent is still —1. However more generally the effort associated 
with approximating Joy- dominates the effort associated with the other two integrals. 

8. Efficient quadrature bases. When multiple Stratonovich integrals contain 
repeated indices, are they as cheap to compute as the corresponding lower dimensional 
integrals with an equal number of distinct indices (none of them repeated)? 

Let i ■ ■ - ip denote the multi-index with p copies of the index i. Repeated integra- 
tion by parts yields the formulae 

9 r 

Ji---ipji---iq = ^ ' (~-i) Ji---ikJi---ip3i---iq-k +(^1)''^ / Ji-ipJi-iq'^Jj i (8.1a) 
k=l •' 
9-1 ^ 

J^-^.J-Jq = YS-^^^-^r-i^^^-^^r-H-. + (-1)'+' / •^-.dJ,...,, . (8.1b) 

k=l 

The first relation (j8.1ap suggests that any integral of the form Ji-.-i^ji-.-i^ can always be 
approximated by a single sum. This last statement is true for q — I. If we assume it 
is true for q—l and apply the relation (|8.1ap we establish by induction that Ji...i 
can be approximated by a single sum. A similar induction argument using (|8.1bp then 
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also establishes that any integral of the form Ji---ipj---jg can also be approximated by 
a single sum. Hence in both cases the quadrature effort is proportional to QN. 

Implicit in the relations (|8.ip is the natural underlying shuffle algebra created 
by integration by parts (see Gaines [20l [19], Kawksi [32] and Munthe-Kaas and 
Wright [48]). Two further results are of interest. Firstly we remark that by inte- 
gration by parts we have the following two shuffle product results: 

If we replace {ii, 12, «3, 14} by j} in (|8.2bp and (|8.2ap and then by 

in (|8.2ap and (|8.2b[) . respectively, we obtain the linear system of equations 



/I 


1 


1 










^JiiJjj Jiijj Jjjii\ 





2 





1 




Jijji 
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1 




Jjiij 




JiijJj '^'^iijj 


Vo 








V 




\Jijij / 




\ Jij Jij — ^ Jiijj / 



(8.3) 



By direct inspection the coefficient matrix on the left-hand side has rank 4 and so 
all the multiple Stratonovich integrals Jjiji, Jijji, Jjnj and Jijij can be expressed 
in terms of Jujj and Jjjn and products of lower order integrals, all of which can be 
approximated by single sums (note that Jiji = JijJi — 'ZJnj^. 
Now consider the set of multiple Stratonovich integrals 

iJ = {Jiii2i3i4is '■ { 

where we exclude the elements Jmjj and Jjjm which we know can be approximated by 
single sums from (|8.1b[) . By considering the shuffle relations generated by products of 
the form. Ji^i2i3i4,Ji^, "^21*2*3 "^^4^51 '-^^1*2 '^^3^4^5 "^^s and Ji^Ji^Ji^Ji^^Ji^ and substituting 
in the 10 elements with indices from 'permsji, i, «,_;/, j}' we obtain an linear system 
of equations analogous to (|8.3p with 50 equations for the 8 unknowns in Z- However 
direct calculation shows that the corresponding coefficient matrix has rank 7. In 
particular, all of the multiple integrals in -3 can be expressed in terms of Jmjj , Jjjm 
and Jjijii- Hence the set of multiple integrals with indices from 'permsji, j, i, j, j}' 
cannot all be approximated by single sums, but in fact require a double sum to 
approximate Jjiju. 

For simplicity assume d — I. Consider numerical schemes of increasing order M. 
If 3/2 < M < 3 all the necessary multiple integrals can be approximated by single 
sums — at the highest order in this range indices involving permutations of {1, 1, 1, 1, 0} 
and {1, 1, 0, 0} are included for which the corresponding integrals can be approximated 
by single sums. For methods of order M > 7/2 we require at least double sums to 
approximate the necessary multiple integrals. 

When d = 2, for methods of order M ^ 1,3/2 the integrals involved can be 
approximated by single sums, but for M = 2 integrals involving indices with permu- 
tations of {2, 1, 0} are included which can only be approximated by double sums. If 
there were no drift vector field then for 1 < M < 2 the necessary multiple integrals 
can all be approximated by single sums, but for M — 5/2 we need to include multiple 
integrals involving indices with permutations of {2, 2, 1, 1, 1} which require approxi- 
mation by double sums. We can in principle extend these results to higher values M, 
however methods of order M > d/2 for d > 3 are not commonly implemented! 

9. Numerical simulations. 



(8.2a) 
(8.2b) 
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9.1. Riccati system. Our first application is for stochastic Riccati differential 
systems — some classes of which can be reformulated as linear systems (see Freiling [TH] 
and SchifF and Shnider [50J ) . Such systems arise in stochastic linear-quadratic optimal 
control problems, for example, mean-variance hedging in finance (see Bobrovnytska 
and Schweizer [5] and Kohlmann and Tang [34]) — though often these are backward 
problems (which we intend to investigate in a separate study). Consider for example 
Riccati equations of the form 



^ d rt 

Ut — Uq + 

i=0 ' 



J2 I {UrMT)Ur+B^{T)Ur+UrC^{T)+D,{T))dWl. 

li y = {U V)'^ satisfies the linear stochastic differential system (jl.l[) . with 

/B.(t) D,it)\ 

then u = UV~^ solves the Riccati equation above. 

We consider here a Riccati problem with two additive Wiener processes, and 
, and coefficient matrices 

D,^(l f), D,^(' U and D,^(\ 0, (9.1) 



2 200/ \ 2 



and 



^0 = ( } \ ) , and Co = ' " 



2 



2 



All other coefficient matrices are zero. The initial data is the 2x2 identity matrix, 
i.e. Mo = I2 and therefore Uq = I2 and Vq — I2 also. We found Higham [28] a very 
useful starting point for our Matlab simulations. 

Note that for this example the coefficient matrices ai and 02 are upper right block 
triangular and therefore nilpotent of degree 2, and also that 0102 and 0201 are iden- 
tically zero so that in particular [ai, a2\ = 0. The number of terms in each integrator 
at either order 1 or 3/2 is roughly equal, and so for a given stepsize the uniformly 
accurate Magnus integrators should be more expensive to compute due to the cost 
of computing the 4x4 matrix exponential — we used a (6, 6) Fade approximation 
with scaling to compute the matrix exponential. See Moler and Van Loan ^47j and 
also Iserles and Zanna [31] , the computational cost is roughly 6 times the system size 
cubed. Also note the order 1 integrators do not involve quadrature effort whilst the 
order 3/2 integrators involve the quadrature effort associated with approximating Jio 
and J20- For comparison, we use a nonlinear Runge-Kutta type order 3/2 scheme for 
the case of two additive noise terms (from Kloeden and Platen [331 P- 383]) applied 
directly to the original Riccati equation: 

^t„,t„+i = 5t„ + f{st,:)h+ DiJi + D2J2 

+ + fiyn + f{y^) + ./(v) - 4./(^t„)) 

+ i7h ((^(^1^) - /(^r))-^io + (/(^2+) - i{y2))J2^) , (9.2) 



where ^ St„ + |/(5t„) ± D^Vh and f{S) = SAqS + BqS + SCq + Dq. 
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Fig. 9.1. Global error vs stepsize (top) and vs CPU clocktime (bottom) for the Riccati problem 
at time t = 1. The Magnus integrators of order 1 and 3/2 shown are the uniformly accurate Magnus 
integrators from Section\4\ 

In Figure 19.11 we show how the global error scales with stepsize and also CPU 
clocktime for this Riccati problem. Note that as anticipated, for the same step size 
(compare respective plot points starting from the left) , the order 1 Magnus integrator 
is more expensive to compute and more accurate than the order 1 Neumann integrator. 
Now compare the order 3/2 integrators. For the nonlinear scheme (|9.2p . we must 
evaluate f{S) five times per step per path costing 20p'^ + 54p^ flops — here p — 2 refers 
to the size of the original Riccati system. For the Neumann and Magnus integrators 
the evaluation costs are 16(2p x p) = 'S2p'^ and 6(2^)^ + ll(2p)2 = 48p^ + 44^^ flops, 
respectively (directly counting from the schemes). Hence for large stepsize we expect 
the Neumann integrator to be cheapest and the Magnus and nonlinear Runge-Kutta 
integrators to be more expensive. However for much smaller stepsizes the quadrature 
effort should start to dominate. The efforts of all the order 3/2 integrators will not 
be much different and the Magnus integrator then outperforms the other two due to 
its superior accuracy. 

9.2. Linear system. Our second numerical example is for a homogeneous and 
constant coefficient linear problem involving two Wiener processes with coefficient 
matrices = Di where the Di, i = 0,1,2 are given in (|9.ip and do not commute, and 
initial data yo = (5 1)"^- In Figure [9721 we show how the error scales with stepsize 
and CPU clocktime. We see that the superior accuracy of the Magnus integrators 
is achieved for the same computational cost. Note that in the case M = 1 we have 
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Fig. 9.2. Global error vs stepsize (top) and vs CPU clocktime (bottom) for the model problem 
at time t = 1 with 2 driving Wiener processes. The error corresponding to the largest step size takes 
the shortest time to compute. 




Fig. 9.3. Confidence intervals for the global errors of the uniformly accurate Magnus integrators 
for the model problem at time t = 1 viith 2 driving Wiener processes. 



d = £ = 2. For the case when h < her, when computational cost is dominated by 
quadrature effort, the relation between the global error £i and computational cost U 
is, ignoring Z//°™' and taking the logarithm. 



\og£i « logXi - ^\ogU 
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Fig. 9.4. Global error vs stepsize (top) and vs CPU clocktime (bottom) for the model problem 
at time t = I with 3 driving Wiener processes. The error corresponding to the largest step size takes 
the shortest time to compute. 

For the order 1/2 Magnus method the computational cost U is given solely by the 
evaluation effort and therefore we have 

log f 1/2 = log/4ri/2 + ^\ogT{ci/2p'^ + ce) - ^logZ^. 

Further note that and Ki are strictly order one, and that T — 1, p ~ 2, ci/2 — 5 
and ci = 7 (counting directly from the inegration scheme). In addition ce ^ Gp'^ = 48 
flops using the (6, 6) Fade approximation with scaling. Substituting these values into 
the difference of these last two estimates reveals that 

logfi -logfi/2 « -i log 68 » -0.9, 

which is in good agreement with the difference shown in Figure 19.21 

Also for this example, we have shown 90% confidence intervals in Figure [?751 for 
the global errors of the uniformly accurate Magnus integrators of orders 1 and 3/2. 
The confidence intervals become narrower as the stepsize decreases and order of the 
method increases, as expected — see Kloeden and Platen [33l pp. 312-316]. 

In Figure 19.41 we consider a linear stochastic differential system driven by three 
independent scalar Wiener processes, with the same vector fields as in the linear sys- 
tem with two Wiener processes just considered but with an additional linear diffusion 
vector field characterized by the coefficient matrix 
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which does not commute with ag, ai or a2. Using Table 7.1 we would expect to see 
for small stepsizes for the order 1 and 3/2 Magnus and Neumann methods, that the 
global error scales with the computational effort with exponent —1/2. This can be 
seen in Figure [9^ 

10. Concluding remarks. Our results can be extended to nonlinear stochastic 
differential equations. Consider a stochastic differential system governed by {d + 1) 
nonlinear autonomous vector fields Vi{y) — instead of the linear vector fields ^ai{t)y^ 
in p.ip — and driven by a d-dimensional Wiener process (W-^, . . . , W^). If we take the 
logarithm of the stochastic Taylor series for the flow-map we obtain the exponential 
Lie series (see Chen [15] and Strichartz [54] ) 

d d 
i=0 j>i=0 

Here [• , •] is the Lie-Jacobi bracket on the Lie algebra of vectors fields defined on W. 
The solution yt of the system at time < > is given by yt = exp at ° yo (see for example 
Ben Arous [4] or Castell and Gaines [HI [14]). Across the interval [i„, t„+i] let at^,t^+i 
denote the exponential Lie series truncated to a given order, with multiple integrals 
approximated by their expectations conditioned on intervening sample points. Then 

yt„+i = exp(af„,f„+i) o yt^ 

is an approximation to the exact solution yt„^-i at the end point of the interval. The 
truncated and conditioned exponential Lie series o't„,t„+i is itself an ordinary vector 
field. Hence to compute yt„^i we solve the ordinary differential system 

u'{t) = (Tt„^t„+i O u(t) 

with m(0) = yt„ across the interval r e [0, 1] (see Castell and Gaines [13l [14]). If 
we use a sufficiently accurate ordinary differential integrator commensurate with the 
order of the truncation of the exponential Lie series then u(l) w yt„+i to that order. 

Hence the order 1/2 exponential Lie series integrator is more accurate than the 
Euler-Maruyama method. Further in the case of commuting diffusion vector fields, 
the uniformly accurate exponential Lie series integrators of order 1 and 3 /2 are more 
accurate than the stochastic Taylor approximations of the corresponding order. This 
generalization is discussed in Malham and Wiese [42j . An important future investiga- 
tion to justify the viability of these schemes would be the relation between global error 
and computational cost, which must take into account the additional computational 
effort associated with the ordinary differential solver. 

An important application for our results that we have in mind for the future are 
large order problems driven by a large number of Wiener processes. High-dimensional 
problems occur in many financial applications, for example in portfolio optimization or 
risk management and in the context of option pricing when high-dimensional models 
are used, for example for the pricing of interest rate options or rainbow options. 
Large order problems also arise when numerically solving stochastic parabolic partial 
differential equation driven by a multiplicative noise term which is white noise in 
time and spatially smooth (see for example Lord and Shardlow [37]). Here we think 
of projecting the system onto a finite spatial basis set which results in a large system 
of coupled ordinary stochastic differential equations each driven by a multiplicative 
noise term. The high dimension d of the driving Wiener process will now be an 



24 



Lord, Malham and Wiese 



important factor in the computational cost for order Af > 1 as for example we will 
need to simulate ^d{d— 1) multiple integrals J^-; though the results of Wiktorsson [5S] 
suggest this can be improved upon. Krylov subspace methods for computing large 
matrix exponentials would be important for efficient implementation of our methods 
for this case (see Moler and Van Loan [47] and Sidje [52]). 

Lastly, extensions of our work that we also intend to investigate further are: (1) 
implementing a variable step scheme following Gaines and Lyons [22], Lamba, Mat- 
tingly and Stuart [35], Burrage and Burrage [TU] and Burrage, Burrage and Tian [S] — 
by using analytic expressions for the local truncation errors (see Aparicio, Malham 
and Oliver [l]); (2) pricing path-dependent options; (3) deriving strong symplectic 
numerical methods (see Milstein, Repin and Tretyakov [15]); and (4) constructing nu- 
merical methods based on the nonlinear Magnus expansions of Casas and Iserles 
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Appendix A. Neumann and Magnus integrators. We present Neumann 
and Magnus integrators up to global order 3/2 in the case of a d-dimensional Wiener 
process {W^, . . . , W^), and with constant coefhcient matrices qq and Ui, i = 1, . . . ,d. 
The Neumann expansion for the solution of the stochastic differential equation (|l.ip 
over an interval [i„,t„+i], where i„ — nh, is 

2/w„+i « + Si,2 + S, + 53/2) 2/t„ , 
where (the indices i, j, k run over the values 1, . . . , d) 

Si/2 = Joao + ^ Jtai + ^ Jiittf , 

i i 

Si — ^ ^ Jijdidj , 
'S'3/2 = ^("/ioooai + Jaidiaa) + ^ JkjiaiUjak 

+ ^^(^OMfliflo + Jiioo-oO'i) + Jajjofoj + JoaaQ . 

i i,j 

The corresponding approximation using the Magnus expansion is 

mag t , , \ 

Vt^L+i ~ exp(si/2 + si + 53/2) yt„ , 
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where, with [•, •] as the matrix commutator, 

i 

Si=Y^ \{Jji - Jy)[a,,aj] , 

i<j 

S3/2 = ^("^jO — Joi)[o-0, O-i] + ^^(Jiij — \ JiJij + tl '^J^ ) ['^i ; ['^i ; %]] 

+ ^ {{Jijk + \JjJki + \JkJij — "^JiJjJkjlaijlaj^ak]] 

i<j<k 

+ {Jjik + \ JiJkj + \ JkJji — ^JiJjJk)[aj, [ai,ak]]) 

i 

To obtain a numerical scheme of global order M using the Neumann or Magnus 
expansion, we must use all the terms up to and including Sm or sm, respectively. 
Leading order terms of order M + 1/2 with non-zero expectation can be replaced 
by their expectations (as detailed at the end of Section [3]). Further, extending these 
solution series to the non- homogeneous and/or the non-constant coefficient case is 
straightforward. 
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